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ORBITAL MOTION OF THE SOLAR POWER SATELLITE 


by 

Otis F. Graf, Jr. 


1.0 Introduction 

It has been proposed to put a series of large satellites 
into geosynchronous orbit for the purpose of collecting solar 
energy and redirecting it toward the earth via microwave radi- 
ation. Preliminary studies are being carried out at JSC on 
the feasibility of these Solar Power Satellites (SPS). 

The large area of the collecting surface (approximately 
144 square kilometers) means that solar radiation pressure 
will cause significant perturbations on the SPS orbit. In 
fact solar pressure will be as important as gravitational 
perturbations. This report documents a study on the effects 
of solar radiation pressure on the SPS orbit. It will be 
shown that the eccent' tv of the orbit can get rather large 

(.08) even though i . litially zero. This is the primary 

difference between oPS orbi*" and other geo.synchronous 

satellite orbits. 

The SPS configuration being considered here is described 
in a study report by the Johnson Space Center (Reference 1). 
Others are discussed in References 2. 3 and elsewhere. How- 
ever, the results in this report are applicable to any geo- 
synchronjus satellite that resembles a flat surface that con- 
tinually faces the sun. 

The main purpose of this report is to investigate the 
orbital evolution of the SPS over its expected thirty year 
lifetime. As a first step, it is assumed that the satellite 
is in free flight, i.e. 'ere is no active orbit control. 
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This will make evident the important orbital motions. One 
of the goals of this study is to describe the motion with 
analytical formulas. These could then be used as a basis 
for developing an orbit control theory that will minimize 
station Keeping costs. 

The perturbing forces cting on the satellite are dis- 
cussed in the next section. To a first approximation, three 
types of forces can be considered separately since they have 
different effects on the orbit. 

(1) Longitude dependent tesseral terms in the earth’s 
geopotential cause a slow drift of the satellite's 
mean longitude. 

(2) Sun and moon gravity cause a rotation of the 
orbital plane. 

(3) Solar radiation pressure will cause an increase 
in the orbital eccentricity. 

Variations in orbital eccentricity e are discussed in 
Section 3. Analytical solution methods are used to develop 
equations for the variation in eccentricity and argument of 
perigee as a function of time. These equations are valid for 
arbitrary initial values of eccentricity and inclination. It 
is shown that e will have a periodic variation with an am- 
plitude of .04 and period of one year. There is also a linear 
increase so that e will grow to .08 within thirty years. 

Earth-Sun-Moon gravity will cause long period variations 
in e . These effects have been studied with numerical inte- 
gration methods and are discussed in Section 4. Evolution of 
the orbital elements is shown for a variety of initial con- 
ditions. The maximum value of e can be reduced by an ap- 
propriate choice of initial conditions. 

Implications of non-circular, non-equatorial geosynchro- 
nous orbits for the SPS are discussed in Section 5. It is 
shown that the daily wariation in longitude is 2e radians. 
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However, these orbits offer certain advantages for the SPS 
and should be further evaluated for their impact on the energy 
collection, transmitting and receiving systems. 

2.0 Perturbing Forces 

The perturbations due to sun-moon gravity and non-spher- 
ical earth have been extensively discussed in the literature 
and only an overview will be given here. .Aicceleration due to 
solar radiation pressure will be derived in this section, con- 
sidering the expected physical dimensions of the SPS. 


2.1 Non-sphericity of the Earth 

This perturbation arises from the fact that the earth is 
not symmetrical about its spin axis. A slice of the earth 
perpendicular to its spin axis has an almost elliptical shape. 
Since the earth rotates once a day and the satellite makes 
one revolution in approximately one day, these gravitational 
perturbations act in the same direction over a long period 
of time. As a result there is a large, long period drift in 
the geographic mean longitude of the satellite (Reference 4). 
The other orbital elements are not sc severly affected. Ref- 
erences 5 and 6 give a good description of this motion. 


2.2 Luni-Solar Gravity 

The luni-solar perturbations have a substantial effect 

upon the node h and inclination I of the orbit. Coupling 

between the .sun, moon and earth’s oblateness (^ 2 ^ cau;e 

large, long period perturbations in I (Reference 7). Table 

I shows some representative values of the inclination after 

two and 26.5 years. If = 0 . the inclination grows to 

14.7® after 26.5 years. An important case is when = 7.3® 

and h =0. Then the inclination and node are almost constant, 
o 
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TABLE I.- VARIATION OF INCLINATION 


I 

0 

h 

0 

I (2 yrs.) 

Max. I . (26.5 yrs. i 

n 

undef . 

1.73° 

14.7° 


270° 

.74° 

14.9° 

B 

o 

o 

0^ 

o 

o 

o 

CM 

15.0° 

mm 

180° 

o 

o 

o 

00 

29.4° 

B 

0° 

7.30° 

1 

7.3° (const.) 


2.3 Solar Radiation Pressure 

The magnitude of the solar radiation pressure depends on 
the weight and cross sectional area of the satellite. The 
most important effect is a rotation of the line of apsides 
and a periodic variation in the eccentricity with a period 
of about one year. To compute the perturbing acceleration, 
the following assumptions are made; 

(1) The SPS is a flat plate cf 10'?<' reflectivity. 
(Reference 1, Section TV.B.l). 

(2) The flat plate maintains an in(?rtial orientation 
perpendicular to the satel 1 i te-sun line. 

(3) Pressure from the microwave transmission can be 
neglected. 

(4) The eartir.s shadow can bo ni‘glecu>d. 

The solar radiation on a flat plate in the v'icinity of 
the earth is (Reference 8): 

9.02-10 N/m' ( lOOv refl<H*ting body) 


•l.51-10~ N/m^ (Blackbody) 
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Thus , for a lO'o reflecting body, the solar pressure is 

4. 96- 10"® N/m^ . 

Solar array area and weight ranges are given in Reference 
1, Figures IV. A. 5. 2, 

97 km^ < area < 186 km^ , 

48 • 10® kg < weight < 123 . 10®kg . 

For the analyses carried out in this report, the following 
"nominal" values were taken: 

Area = 143 km^ 


Weight = 82.5 • 10‘' kg. 

Let F be the force due to solar radiation pressure and 
M the spacecraft weight. Then the perturbing acceleration is 



M 


where the magnitude A=1 a| is constant. Let S be the sur- 
face area in square meters, then 




( 2 . 1 ) 


where <j is t!ie acceleration of gravity the surface of 
the earth (;7=9.807 m(sec)”' ). If M is expressed in kilo- 
grams, then 


- 5 .06 • 10 


S 

M 


( 2 . 2 ) 
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A 

Note that - Is unitless. Using the nominal values for 
9 

area and weight, 


and 

1 » 1.73 m2 kg"' , 

(2.3) 


- = .875 • 10"® 
9 

(2.4) 

Taking 

into account the expected * nge in size 

A 

and weight, - 

can be 

in the interval 


. 72 • 10" ® < - < 1 . 16 • 10" ® . 

9 

(2.5) 


One additional comment needs to be made on assumption 
(4). An equatorial geosynchronous satellite will pass 
through the earth's shadow once a day during the eleven days 
before and after the equinoxes. It will re. iain in the shadow 
for a maximum of 75 minutes on the day of th quinox. Th“ 
amount of time in one year that the satellite is in the shad- 
ow is small and will not be important in .studying the long 
term effects of solar radiation pressure. 


3.0 Solar Radiat io n Pressure Effects on the O r bit 

Variations in orbital eccentricity due to the perturbing 
effects of solar radiation pressure are discussed in this sec- 
tion. The magnitude of the perturbing acceleration was dis- 
cussed in Section 2.3. An approvimate solution is given for 
the variation of e as a function of lime. This .'solution is 
valid for small eccentricities, i.e. e .08 Comparison 

to numerical integratiiin shows that the solution is valid for 
about eight years. After that time. gravitativ)nal effects 
(discussed in Section 4) become important nowever. this so- 
lution shows the geiieral nature of the p.-=‘r*'urbations in 
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eccentricity and argument of perigee. Also, it could be use- 
ful to compute station keeping maneuvers for orbit control 
purposes . 

Musen (Reference 9) did some early work on orbit per- 
turbations due to solar radiation pressure. He was concerned 
with the orbit of Vanguard I where the rotation of the line 
of apsides (due to oblateness of the earth, J 2 ) was nearly 
conmensurate with the motion of this sun. This caused large 
perturbations in the height of perigee. Hori (Reference 10) 
developed a canonical theory for this resonance problem. 

Solar pressure was assumed by Musen and Hori to be order of 
magnitude * 

The case where solar radiation pressure is large (such 
as with the SPS) has been treated by Zee (Reference 11), 

Bosch (Reference 12), Ahmad and Stuiver (Reference 13), and 
Van der Ha and Modi (Reference 14). The analyses of Bosch 
and Ahmad and Stuiver are restricted to motion in the ecliptic 
plane with the sun assumed fixed. Their results are thus 
valid for only a few revolutions of the satellite. Zee shows 
that the eccentricity will have a period of one year, but he 
considers only the case where e is initially zero, and 
does not give any quantative results. Van der Ha and Modi 
use the two variable expansion procedure to describe the 
yearly motion of e for the case where the orbit lies in the 
ecliptic plane. They use an area to weight ratio of 20, where- 
as References 1, 2 and 3 indicate a value near 2 or 3 (see 
equation 2.3). These investigators did not consider the im- 
portant secular increase in eccentricity or coupling between 
radiation pressure and gravitational pertrubations . 


3.1 T..e Solar Radiation Perturbing Function 

Let r be the satellite's position vector referenced to 
an earth-centered coordinate system whose x-axis is in the 
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directiOQ of the earth's north pole. The x- and y- axes 
lie, therefore, in the equatorial plane. 

1 * 

The acceleration vector of the satellite is 

** p 

r = A - + ~ . (3.1) 

p 3r 

* ft 

where D is the gravitational force function . p is the 

vector from the sun to the satellite (Figure 1). 


Sun 



Satellite 


FIGURE 1: Earth-Satellite-Sun Geometry 


f 

Dots refer to derivatives with respect to time, i.e. 



++The sun, moon and earth gravitational effects ere included 
in U* . 
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Let be the vector from the earth to the sun. De- 


fine the unit vector 


6 ■ I * - i <^<9 - ?) 


p • lpl 


r@ « and r » |r| 


For a geosynchronous satellite in a nearly circular orbit. 


r - 42»164 km. 


The earth-sun distance is 


Tq - 149. 5 *10 km. 


Therefore, the ratio — will be small, i.e, 


Ir = 2.810"*’ 


The unit vector & can be expressed in powers of the 
small parameter. From the law of cosines (see Figure 1): 


p* = r^ + r - 2 r^ r i|i , 


1 = — fl + (— 

P L 


) - 2( — ) aos <{» 


The above expression can be expanded in powers of — 


- = ^ Z (— )" P (c?os <!/) 
P o r^. n 



where P^(C) is the Legendre polynomial with argxuaent C • 
The expression for 0 is then 

Cl 5 

From the above expression it is seen that the replacement 


■* 

P 

P 



(3.5) 


involves an error of 2. 8 •10**'* . The equations of motion 
are then 



(3.6) 


The components of (3.6) in rectangular coordinates are 



X • - A — 



+ 


•■o 

dx 



-♦ 



9U 


y = - A -2. 

+ 


*0 

3y 



au* 


z - - A 




dz 

Define the 

new force function 

II 

*r !© 

- A x — + y -S. ^ 






(3.7) 


(3.8) 


The differential equations of motion are then 
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For perturbation problems, it is desirable to write the force 
xunctioD in the form 


U - ^ (1-V) , (3.10) 

where 

y - y, Vq* \ + V, (3.11) 

is the "perturbing function", u is the gravitational con- 
stant for the earth (3.98601-10® km® sec is the 

contribution of solar radiation pressure and can be written 
as 



A similar perturbing function was used by Hori (Reference 10). 

represents the geopotential. and Vq are the grav- 

itational potential functions of the moon and sun, respectively 
(Reference 7). 


3 2 Order of Magnitude Considerations 


This section considers the magnitudes of the 
terms in (3.11). It is shown in Reference 7 that 
synchronous satellite in a nearly circular orbit, 
tudes of the gravitational terms are 


various 
for a geo- 
the magni- 


|v„| = 2.5-10' 

|Vj I = 1.6-10 


|Vg| = .75-10 
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From equation (3.12), 



(3.13) 


where a is the semi -major axis. 


It > as shown in Section 2.3 that a typical value of 


is 

Since 


A - .875*10"® g 


V 



A 


where is the radius of the earth, the magnitude of 

can be written as 


or 



A 

ff 


a 

Rw 


» 


|Vg| - 3.84*10"® 


Therefore, solar radiation pressure has the order of magnitude 
of the gravitational terms. ' 

Since V^, V^, and have nearly the same mag- 

nitudes, it is allowed, for a first approximation, to consider 
each effect separately in arriving at an analytical solution. 

The following section will, therefore, investigate the pertur- 
bations in the orbital elements due to solar radiation pressure. 


3.3 Delaunay-Similar Elements 

The solution will be developed through the use of canon- 
ical element differential equations in an extended phase space. 
Delaunay-Similar elements in the eccentric anomaly (DSu) have 
been presented in References 15 and 16. The angular variables 


are: 
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u » eccentric anomaly, 
g “ argument of perigee, 
h « argument of the ascending node, 
t » time element. 

The action variables are: 

U * related to the two-body energy, 

G = total angular momentum magnitude, 

H = z-component of the angular momentum, 
L = negative of the total energy. 

4* 

Differential equations for these variables are 


da^ 3F 
dx 3Bj 


d6. 3F 

— = . (i=l,2,3,4) (3-14) 

dt 30j^ 


where the Hamiltonian function is 


F = U ^ V (3.15) 

/2C 


and V is the perturbing potential function. The time is 
given in terms of DSu-elements by the equation 


t 


t 


2L 


e sin 


u 


(3.16) 


In unperturbed motion. 


u 

Z 


T + constant 


P 

(2L)^* 


T + constant 


(3.17) 


t 

The following element notation is used: 


a 

1 

= U, 

a 

2 

= g. 

a = Ij , 
3 

a - Z , 

4 

6 = II. 

1 

B = 
2 

6 

= H, 

B 

« L. 






3 4 



- 22 - 


and the remaining elements are constants. 

The following abbreviations are used (Reference 16): 



U* 


✓2D 


£ B 1 . e 008 u , ein 1 • /. 

mr 

(3.18) 

a 

G* 

is numerically equal to zero. 

so that 

U is defined 

U • [l - ^ 

• 

(3.19) 

✓2D L -1 



Therefore, since both U and L depend on V,, e and a 
are slightly different from the Instantaneous eccentricity 
and semi-major axis, respectively, in the case of perturbed 
motion (V*0). 


3.4 Development of F in Terms of Elements 

Considering only perturbations due to solar radiation 
pressure, the hamiltonian is 


where 


F ■ U H-. + p 

✓2D ® 


F - V 
* ^2D ® 




Using (3.12), 



(3.20) 
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The coordinates are given in tenns of elements by 
* - a[f,^(oo« u-e) + Cg ein 

y » aj^ ((?os u-e) + /1-e' ein ^ , 

z » a|n^(<3oa u-e) + /1-e ein , 


where the following abbreviations are used: 


* COB h 008 g - ein h ein g ooe I , 


C - ooe h ein g - ein h ooe g ooe I , 
2 


= ein h ooe g + ooe h ein g ooe I , 


C = sin h ein g + ooe h ooe g ooe I 
2 


n = sin g sin I 
1 


n = cos g ein I , 
2 


The direction cosines of the sun are (see Appendix A): 


— = g^) - e^j cos g^ 


'© 


— = C sin{t +g^) -Ce sin g. 


~ = S sinil^+g^) 


S sin g^ 
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where the following notations are used: 


C - oo«(23®27') 

5 * «i«(23®27') 
n t + 

(D O Oo 

n » 360® ' (365. 2422 days) 
G) 

e^ = .0167 . 

= 281* 

o 


0 


n 

o 


t + A 


eo 




(3.24) 


When expressions for x,y,z,r and equations (3.23) are 
inserted into (3.20), is given in terms of elements: 


Aa* 
F = 


||(l+e^) eo8 u - ^ e - ^ e aoe 2i^ 


• jc^(eo8 9 - <^oe g^) + (atn 9 - ein g^^jj 

r 1 -1 (3.25) 

+ (1-e*) I ein u - •g’ e ein 2u • 


0 


^ (ooe 6 - 6q 008 g^) + N (sin 0 - ein 


& 


8®)J 


where 


N'CC+5n , N-C?+5n 

111 222 


. (3.26) 







The variable 6 in (3.25) contains the time, 


3 ® n t + + gffl 




(3.27) 


Using the time equation (3.16), 6 can be considered an ab- 

breviation involving DSu-elements ; 

« 

n U 

0 “ ^ “ — S— e 8 in u + So • (3.28) 

2L 


Carrying out the products in (3.25) , 


Aa 


F = 


V j(l+e2)C^ - (l-e^) nJ ao8 (0+u) + 


2/2I7J 

+ [^(1+eM + (1-e^) cj sin (0+u) + 

+ l^l+e^) + (1-e^) N^Jeos (0-u) + 


+ (1+e^) N - (1-e^) C sin (0-u) - 
_ • ?J 


- 3 e ? 008 0 - 3e N sin 0 + 

1 1 

♦ ge(l-e“) - |e c] 008 ( 0+2u) - 

e|N + (1-e^) (0+2u) - 

-| e|j;^ + (1-eM nJcos (0-2u) - 

“ (1-e^) (0-2u) - 


(3.29) 


[ 3 1 

(l+e^) 008 u - ^ e - ^ e oos 

- 2 cos 8 in ^)(l-e^ ) |s£n u - ^ e ain 2i^| . 
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Slnce the Interest is in long period motion, short pe- 
riod terms (those periodic in u) will be eliminated from 
. This can be done at the same tim^ that the time equa- 
tion (3.28) is inserted into F . The terms dependent on 

8 

time in (3.29) are 


ain 

008 


(0 + i u) 


1 « 0 ,± 1 ,± 2 . 


or 


^ ^ + S® + iu). i=0,±l,±2. 


where 


Un e 
X S_ 

2L 


f 

The following relations (Reference 17, p.2 -0) are used 

«|>00 

ein(a + 6 sin y) ® E J. (8) ain(j y + u) 


+00 

aoaia + 6 ooa y) “ ^ (S) ooaii y + a) 

— 00 J 


Thus , 

sin{tiQ t + X sin u + + 


+ i u) 


+00 p 

« E J.(X)8inn^t 

•1 S — oo ^ I— 


(3.30) 



The mean of a function f(u) with respect to u 

r2ir 


<«“)>„ • ^ 


f(u) du 


is defined as 


Jj(6) . j-0. 1. 2,--«, 


are the Bessel coefficients. 
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Therefore, using (3.30) , 

+ i u)>^ - J_^(X) einCng, £ + *<^0 ®o)- (3.31) 

Similarly, 


^c/^8(0 + i u)^^j = A + ^©o (3.32) 


Jsing (3.29), (3.31) and (3.32), the elimination of shor^ 
period terms results in 

Aa^ 


<F > = Jf(l-e*) J (A) N - 

® “ / 2 U ! L 1 2 


] 


" 2 ® ^^*^0 ^ ^ *’®o S®) 


- (1-eM J (A) c + 
_ 12 


(3.33) 


+ I e (3J^(A) + J^(A))nJ 8in(n^ ^ ^ ^&c * «®> 


2 ® 8t« g^) 


Note that several terms cancelled because of the identity 


J.(A) = (-1)^ 


The averaged hamiltonian is now 


/ 

\ 





(3.34) 


Since u does not appear explicitly in • U will be 

a constant. The remaining developments will concern (3.34) 
only. Therefore, the notation ^ will no longer be needed. 
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The hamlltonian in terms of elements Is then 

u 

F - U + F„ . (3.35) 

/2E ® 

where F is a function of g, h, U, G, H, L given by the 
8 

right side of (3.33). 

A further simplification of F^ can be made. Consider 
the Bessel coefficients appearing in (3.33). The argument X 
has been defined as 

Un 

X ■ — ^ e 
2L 


Since L is the negative of the total energy. 


L 




(3.36) 


where a is the Instantaneous osculating semi -major axis. 
From (3.19), 






9 


Inserting (3.36) Into the above expression and expanding in 
t 

powers of ; 


U - 



(3.37) 


From (3.36) and (3.37) , 

^ ^ [l * (X| V,)] , (3.38) 

* 1 * 

In the unperturbed case, U is equivalent to the classical 
Delaunay variable L 



Since n 



is the osculating n^an ration, it is seen that 



is approximately the ratio of the mean ration of the satellite 
to mean motion of the sun, i.e. 



nierefov^, 

1X| s (2.8-10“*) e 
or 

|X| < 2.8-10”® 


(3.39) 


for a geosynchronous satellite. Any tenn depending on X 
will therefore be neglected. 


The Bessel coefficients can be expressed as a power 
series in the argument: 


n 


(X) 


x“ r 

1 + 

2”n! L 2®-l (n+l) 2** l-2 (n+l ) (n+2) 



Since all powers of X can be neglected, Ji(X) and JaCX) 

in F can be set to zero and Jq(X) can be set to one. 

8 

Thus, the expression for F becomes 

s 

Fg = - I G e j^e C^(<?os v - aos g^) + 

+ e Ni(a-Jn V - e^ sin g^)l (3.40) 


where 
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eCi 


eNi 


e eoe g ace h - e sin g ein h aoe I 
C 

* S e ein g ein I , 


e 


e aoe g ein h e ein g aoe h aoe I 


«J 


(3.41) 


V - 1 + g<B . 

Tlie small parameter e Is unitless. (B^nember that A 
has units of acceleration), g is the acceleration of gravity 
at the 8urfc.ce of the earth and the radius of the earth. 

Iherefore , 


But fr^ 
and (3.19) , 


A ^ a P 

e * — (^) 

g Gi^ 

G « U A-e* . 


G/2D = M (1-V„) /l-e^ 
8 


so that 


Since e 


e can be expressed as 

e » - (1-V^r‘ (l-e*)"^ . (3.42) 

g 


is already small, V and e‘ can be neglected 

8 


e 



(3.43) 


Therefore, the small parameter e is the order of magnitude 
of |V^| (see section 3.2). 



-31- 


3.5 EleTODt Differential Equations 

Introduce the non-singular elements 


P • 


P » 



G* 

— oosig-t-h) 
U* 


G 

— 008 h 
H 


/ 72 ? 


U" 


«in(g+h) 


Q 


- /l - - «i 


8xn h 


H 


» 


(3.44) 


These elements are defined for zero eccentricity and inclina 
tion. 

The differential equations for p and q are given by 
by the chain rule: 


dp 

3P 

dG 

3p 

dg 

dp dU dp 

dh 

— — s 

— 

+ 

— + 

— + — 

— 

dx 

dG 

dx 

3g 

dx 

dU dT dh 

dT 


with a similar equation for ^ . Make use of the canonical 

dT 

differential equations 


dg 

- 55 

3F 

dG 

3F 

dT 

3G 

dT 

3g 

dh 

3F 

dU 

3F 

dT 

3H 

dT 

3u 


Also, since F no longer depends on u 


ap 

— « 0 
3u 
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The necessary partial derivatives of p and q are: 


9p 

n* 

aq n* 

ai 

eoa(g4h) , 

— 8in(g4>h) , 

9G 

Ge 

dG Ge 

ap 


dq 

— ■ 

- e 9in(,S*h) , 

— * e eoe(g‘<'h) , (3.45) 

as 


dS 

ap 


aq 

— ■ 

- e 8in(g4‘h) , 

— ■ e eae(g4'h) 

ah 

dh 


The differential equations for p and q' are then 



It is necessary to develop the right sides of equations 
(3.46) in terms of p,q,P,Q . This will be done in the 
following steps: 

(1) Conq>ute the derivatives of F with respect 
to g> G, H . 

(2) Insert these derivatives into the right sides of 
equations (3.46). 

(3) Use equations (3.44) to express the right sides 
of (3.46) in terms of p,q,P,Q . 

These three steps have been carried out in Appendix B. All 
the necessary partial derivatives are given there. The re- 
sulting nonsingular equations are: 



-33 


dp 

dx 




P Q ♦ q Q (qP+Qp) ( COB V • e, 


>]' 


008 g^) - 


-|n*c (i-P*) + n* 5 P (2-(P*+<3*))* - 

- C q P (qP+Qp) ♦ 

+ 5 q (qP+Qi>)(l-/Q*+P*)) (2-(P*+Q*))"^ («i« v - e^sin g^: 

( 3 . 4 *'' 

dq 3 f r -1 

— • — e I |n* (1-Q*) - P Q (qP+Qp^(eo« v - aos ^ + 

-jn* C P Q - n* 5 Q (2-(P*+<J*))* + 

+ C p P (qP+Qp) - 

- S p (qP+Qp) <1-(Q^+P*)) (2-(P*-KJ*))"^(«t« v - e^eos 


Notice that n is a function of p and q 

n = /l-(p*+q*) 

The differential equations for P and Q are developed 
in a similar manner as that done above. The details are car- 
ried out in Appendix B and the equations are shown on the 
following page: 
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dP 3 


* (r 

- - e Q P 
4 L 


P (qO-pP) + 2 pliaoa v - gL ) ♦ 


+ (7p|p (qQ-pP) + 2 i^(«<n v - e^ain g^) + 

+ 5 |l - (P*+Q*)J|2 - (P*+Q*)] *[P (pP-qQ) - 

- 2 ]^(8in V - ain g^) , 

4 

< 

— e QjQ (pP-qQ) + 2 ^(o 0 « v - ooa 6^) + 
+ £?P (pP-qQ) + 2 ^(ai« v - ain g;^) ♦ 

+ S |l - (P='+Q*)J[2 - (P*+Q*)J *[q (pP-qQ) + 
+ 2 ^ (ain V - ain S^)| 


The remaining differential equations come from the ca- 
nonical equations 

dt dF dL aF 

dT 3L ’ dT at 

These are computed in a straightforward manner from (3.35) » 
(3.40) and (3.41): 
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dT 


u 

(2L)V* 


1 * 


9 

- e 
2 




(p*+q*) 


>]Ne 


P (1-Q*) - 


dL 

dT 


- q Q l^(ooa v - ooa + 

+ [c (q-qP*“PPQ) + 

+ S (qP+Qp)(2-(P*+Q*))^(atn v - ooa g^)l 



+ S (qP+Qp)(2-(P*+Q^))^ 008 


(3.49) 


3.6 Orbits in the Ecliptic Plane 

t 

For orbits that lie in the ecliptic plane . 

h = 0 , I = € 

Therefore , 

P = /1-cos r , Q = 0 , 

and the differential equations are greatly simplified: 
dp 3 p -I 

— « e 1 - (p^+q^) (ain v - e^ sin g ) . (3.50) 

dT 2 ^ 

is the angle between the equatorial and ecliptic planes. 
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dq 3 r* -n 

— ■ - e 1 - (p*+q*) (008 V - e- ooe g^) , (3.31) 

dT 2 L -i 


dP 

— - 0 , (3.62) 

dT 


dQ 

— - 0 , (3.53) 

dT 


dl 

dT 



p (008 V - 008 g©) + q (ain v - 


(3.54) 


8 ^« 8 ®) 


0 


dL 3 
— « — n 
dT 2 

- q (008 V - ain g^) 


1 - (p^+q^]^ [p 


p (8%n 


e®c<»8 g^) - 

(3.55) 


Comments 

(1) Equations (3.52) and (3.53) indicate that the orbital 
plane ( l.e. inclination) will remain fixed, as expected. 
This is because there are no out -of -plane perturbations 
on the orbit. 

(2) The equations no longer depend on c 

(3) The equations for p and q contain secular terms that 
are proportional to the eccentricity of the sun's orbit. 
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3.7 Solution of the Linearized Equations 

When the eccentricity and inclination of the satellite 
are small, the element differential equations can be simpli- 
fied by neglecting from equations (3.47), (3.48) and (3.49), 
the second and higher degree terms in p,q,P,Q. The result- 
ing equations are: 


dp ^ r /-I "1 

— e C + 5 / 2 P (ein v - e_ ain g^) 

dx 2 «- -1 


dq 3 p 

— = — e \oo8 V - e_ 

dT 2 L ® 


008 


+ S Q (ein v - e^, ain g^^ , 


(3.56) 


(3.57) 


dP 3 

— - - S c p iain V “ ain , 

dr 4 


(3.58) 


dQ 3 

— = - — /S S e q iain v - ain g^) , 

dx 4 


dt 

dx 


9 p 

1 + — e p (ooa V - e^ ooa g<p) + 
2 L 


(2L)y'^ 

+ C q (ain v - ain 




dL 3 p 

dT 2 ® L 


p sin V - C q cos 




(3.59) 


(3.60) 


(3.61) 


These equations are not linear since v> contains t through 
the equation 
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Also, e depends on L. However, noti^ that the perturba- 
tion In total energy (equation (3.61)) will be small because 
of the coefficient n^ e . It is also periodic. It is 
therefore allowed to let L' be a constant. 

** ' ^ / n 

The perturbed part of the time elen»nt equation (3.60) 
will be small because it is proport ional' to the eccentricity. 
Therefore, let 

U 

I <r — t . (a. 62 ) 

(2L)^ 


The derivatives of P and Q are also proportional to 
e . In fact, the effect of solar radiation pressure on the 
orbital plane is negligible when eos^ared to the gravita- 
tional effects. It is shown in Reference V that sun-moon- 
earth gravity cause a motion of the orbital plane that is de- 
scribed by the following expressions: 


o - a ein (wx + 6^ , 


o 

Q » — a 008 (u)T + 9) 
Y 


(3.63) 


where 


.0902 


1.015 


(I) = 5.170«10"^. 


a and ^ are integration constants that depend on the initial 
values of P and Q. The mean values are 


Q 


(3.64) 


and correspond to the equilibrium solution 


I - 7.31' 


h = o 
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For most initial condition, P and Q will never be far away 
from their equilibrium values. Therefore, insert (3.64) 
into equations (3.56) and (3.57), giving approximate equations 
for the derivatives of p and q : 


dp 3 

— e 5 (sin 9 - sin g^) , 

(3.65) 

dq 3 

— ■ - c (eos 0 - 008 g^) , 


where the additional abbreviations have been introduced: 


$ - (7 + So 


6 


0 ■ 6 T + A _ + g^ 
©o © 


(2L)V^ 


(3.66) 

(3.67) 


(2L)^ 

But since is the mean motion of the geosynchronous 

P 

satellite, 

6 = (365.25)"^ . (3.68) 


Equations (3.65) are uncoupled and can therefore be 
immediately solved: 

p * <I> B (oos 0 + T 6 e^ sin ) + C, 

© © 1 


q « ^ (sin 0 - T 6 e^ qo 8 


(3.69) 


where 


3 e 
2 6 


<t> 


(3.70) 
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and C2 are integration constants that depend on the 
initial values of p and q . 

Comments on the Solution 

1. The motion can be represented in a plane with p.q the 
rv)ctangular coordinates . 

2. In all cases, p,q will describe an ellipse whose center 
has a linear translation. The motion around the ellipse 
has a period of one year. 

3. The mean eccentricity will have a linear increase or 
decrease . 

3.8 Numerical Results 

This section will discuss some quantitative and qualita- 
tive results of the solution developed in Section 3.7. First, 
the solution is verified by coinparing it to a numerical inte- 
gration. Then the solution is used to describe the general 
behavior of orbital eccentricity and longitude of perigee. 


3.8.1 Numerical Experiments 

It is necessary to demonstrate that the analytical solu- 
tion and its associated assumptions are valid. This has been 
done by carrying out cumparisons with a numerical solution 
obtained from the STEPR multirevolution program (References 
18 and 19). Since the purpose here is to check out the accu- 
racy of equations (3.69). only solar radiation pressure was 
Included as a force model option in STEPR. The additional 
effects of gravity will be discussed in Section 4. It will 
also be shown there that equations (3.69) give a good approx- 
imation to the complete problem over a period of a few years. 

Comparisons between the analytical solution and STEPR 
are shown in Tables 11(a) and 11(b). The area to weight ratio 
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vas assumed to be 1.73 kg~*(see Section 2.). The Initial 
epoch was noon January 1, 1980. Eccentricity e based on 
the analytical solution was first computed by evaluating p^ 
and from equations (3.60). The algorithm for computing 

p^ and q^ is described in Appendix C. e^ is then obtai.ned 
from 

®A “ • (3.71) 



TABLE 1 1. -ANALYTIC SOLUTION VERSUS NUMERICAL SOLUTION 

% 

(a) ®o “ ® » Sq '^***^®^i**®^ 


Years 

®A 

®N 

% 

% 

% 

Percent 

®A 

Error 

•\> 

®A 

9.6 

.0481 

.0485 

132.7® 

133.7® 

0.8 

0.7 

19.5 

.0610 

.0620 

148.8 

149.6 

1.6 

0.5 

30.1 

.0511 

.0517 

-173.8 

-174.2 

1.2 

0.2 
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TABLE II. - CONTINUED 
(b) • 021 , » -80.6® 


Percent Error 


Years 

®A 


% 

8a 



8a 

9.6 

.0328 

.0337 

153.7® 

155.0® 

2.6 

0.9 

19.5 

.0502 

.0519 

167.4 

168.3 

3.3 

0.5 

30.1 

.0548 

.0553 

-150.9 

-151.9 

0.9 

0.7 

Comments 







1. The 

results 

in Table 

II show 

that the analytical 

solu- 


tlon gives between 2 and 3 digits of accuracy ever a 
period of 30 years. This is sufficiently accurate to 
describe the general behavior of the orbit. 

2. The accuracy of equations ^3.69) subswuntlates the as- 
sumptions that were made in the course of their deriva- 
tion. The same approach can be used to solve the com- 
plete problem (Including gravitational perturbations). 

3. The errors in e^ and g^ do not increase with time. 
Therefore the analytical solution contains all long 
period effects. 

4. Plots of e versus time are shown in Figures 8 and 9, 
respectively. 


3.8.2 Qualitative Description of the Orbit 

The orbital behavior can be described for different 
initial conditions on eccentricity e^ and longitude of 
perigee g^ « Also, the motion will depend on the epoch 
of initialization since the problem depends explicitly on 
time (i e. on the position of the sun in its Imagined orbit 
about the earth) . The example cases considered Lere are 
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shovn in the table below. These cases were chosen so as to 
demonstrate some of the essential features of the motion and 
to illustrate some preferred orbits. 


TABLE III.- INITIAL CCaiOITIONS FOR EXAMPLES 


Case No. 

®o 


Initial Epoch 

1 

.0 

undeffned 

noon. 

1 Jan. 

1980 

2 

.0210 

-80.6® 

noon. 

1 Jan. 

1980 

3 

.0292 

-34.8 

noon. 

1 Jan. 

1980 

4 

.0214 

-67.4 

noon. 

1 Jan* 

1980 

5 

.0 

undefined 

noon. 

3 April 

1980 

6 

.0 

undefined 

noon. 

1 Oct. 

1980 

For each 

case , the 

evolution of 

elements 

p> q 

and e 


is considered over a period of ten years. This is the maxi- 
mum time interval for which the analytical solution is valid 
(see numerical comparisons in Section 4). The area to weight 
ratii was taken to be 1.73 m^ kg~* , as before. Equations 
(3.69) were programmed on the Hewlett-Packard 9810 micro- 
computer, using the algorithm described in Appendix C. The 
HP9810 plotter was used to produce the plots shown in Figures 
2 through 9. 

Two figures illustrate the results for each case: 

(a) p versus q: This shows the motion in p,q-space. 

Also, e and g are the polar coordinates of the point that 
traces out the curve (see equations (3.44)). The direction of 
motion is indicated by arrows. When the curve passes through 
the origin, :he limiting value of g will be the tangent to 
the curve. 
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(b) e versus time: This figure shows the variations 

in e as a function of time in years. 


Discussion of Results 


Case 1 (Figures 2(a) and 2(b)) 

(1) The motion begins from the origin In the p.q- 

plane and in a direction that is approximately 

t 

90^ from the sun's initial n»an longitude 6^ . 

This can be seen from equations (3.69) with 
6 » 1 ; 

e eo6 g = ♦ ^oe (S t + - cos 0^ , 

e sin g (5 T 9^) “ 9^ 

Notice that e = o for t ** o . Then, 


sin (6 T + 6 ) - sin 0 
''' ° 
tan g ® 

aos (5 T + 0^) - 008 


Using I'Hopital's rule, 

g “ “ ctn - tan ( 90* + 9^^ ) . 

For this case, neon on January first corresponds 

to 9^ » - 80.6® . Therefore g^ » 9.4® 

0 o 


The sun's initial mean longitude is 9^»t^ +g^, where g^ 
is measured from the vernal equinox (x-axls) in the eclip- 
tic plane. See also Appendix C. 
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Case 


(2) The center of the ellip^ In the p,q- plane 

is initially at (C^,C 2 ) and moves toward the 
lower left at the rate of per year. 

Notice that this rate depends on the area to 
weight ratio of the satellite by way of the small 
parameter . The direction cosines of the 
motion are (stn g , - aos g^). Using g =282.5®, 

G G G 

the d-rection cosines are (-.9763, -.2166). 

(3) The eccentricity is periodic and returns almost 
to zero after one year. Its maximum in the first 
year is approximately 2^ . The motion of the 
ellipse in the p,q- plane is seen as a nearly 
linear component in the variation of e 

2 (Figures 3(a) and 3(b)) 

% 

(1) The initial values of anc were chosen 

so that constants C, and C.. were both zero. 
The ellipse is centered initially at the origin. 
For a while, the eccentricity is nearly constant 
at the value of (|> . But as the ellipse moves 

away from the origin, the oscillations in e 
increase in amplitude. Just as the origin is no 
longer inside the ellipse, the amplitude will be 

2 • After that, the amplitude stays the 

same, but there is a linear component to the 
variation of e . (See also Figure 9.) 

(2) As lorg as the origin lies inside the ellipse, 

g will circulate once a year. If the origin 
is still inside the ellipse and the path passes 
near the origin, g will change very rapidly. 

When the origin no longer lies inside the ellipse, 
g varies about the mean value of g.^ - 90 ® 

The amplitude of these variations decreases with 
time . 
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Case 3 (Figures 4(a) and 4(b)) 

(1) The initial conditions for this case were chosen 
so that the ellipse would initially pass throu^ 
the origin and its center would move directly 
toward the origin. 

(2) The important result of this choice of initial 
conditions is revealed by the behavior of e as 
shown in Figure 4(b). in fact Figure 4(b) is a 
mirror image of Figure 3(b). 

C&se 4 (Figures 5(a) and 5(b)) 

(1) This case is similar to Case 3 in that the center 
of the ellipse will pass through the origin. 
However, the motion is such that the eccentricity 
is nearly constant for a longer period of time. 

For "cample the maximum value of eccentricity will 
be chan .025 for four and a half years. 

(2) The . tude of perigee g will circulate. 

Case 5 (Figures 3(a) and 6(b) ) 

This case and the next one show the effect of the 
epoch of initialization on the subsequent motion. 

f\, 

e and g are the same as in Case 1, but the 
o o 

epoch is three months later. The ellipse is mov- 
ing directly away from the origin and the linear 
growth component of e is more severe. This case 
is not desirable when large eccentricities are 
to be avoided. 

Case 6 (Figur s 7(a) and 7(b) ) 

Again, the eccentricity is initially zero, but the 
epoch is six months later than Case 5 and nine 
months later than Case 1. However, the motion is 
very similar to Case 3. The center of the ellipse 
in the p,q-plane moves directly toward the origin 










FIGURE 3: Variations in Eccentricity, Case 2 

(a) p versus q 
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FIGURE 5: Concluded 

(b) e versus Time (years) 
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F GURB 7: Variations in Ec r en t r 1 c 1 1 y , Case 6 

(a) p versus q 





Time (years ) 

FIGURE 9: Long Term Variation in Ec cen t r 1 c 1 ty * Case 2 
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and the amplitude of oscillations in e decreases 
to zero after ten years. 

Figures 8 and 9 show ...e long term variations in e 
for Cases 1 and 2, respectively. Of interest here are the 
effects of the sun's orbital eccentricity e^ over 30 
years. In fact, this secular Increase becomes the dominant 
effect on the motion. However, the numerical studies that 
are discussed in the next section show that gravitational 
perturbations eventually become significant and can actual- 
ly limit the secular growth in e due to e^. 


4.0 Long Period Variations in Eccentricity and Inclination 

The long period changes in the shape and orientation 
of the SPS geosynchronous orbit are discussed in this sec- 
tion. Gravitational and solar radiation perturbations are 
included in the analysis. An analytical solution does not 
yet exist when the problem contains all perturbations simul- 
taneously. Therefore, the results discussed in this section 
are based on a numerical integration of the satellite equa- 
tion of motion. 

As indicated in Sections 1 and 2, the perturbing 
effects can generally be separated as follows: 

(1) Rotation of the orbital plane 

The combined effects of sun and moon gravity and 
the oblateness of the earth cause large, long 
period changes in the inclination. If ini;,ially 
zero, the inclination will increase at the rate 
of .859 degrees per year. Solar radiation pres- 
sure has a very small effect on the inclination. 

(2) Variation in the orbital ea aentriai t'j 

The eccentricity can have large changes, primari- 
ly due to solar radiation pressure. However, 
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gravitational perturbations due to luni-solar 
gravity and earth's oblateness can be iBg>ortant 
over a long period of time.' 

(3) Drift of tho eateltite *8 mean longitude 

This effect is caused primarily because the sat- 
ellite's orbital revolutions are in resonance 
with the daily rotation of the earth. Luni-solar 
gravity has a small direct effect on the drift in 
longitude. There is an indirect influence due to 
solar radiation pressure ( e varies ) and luni- 
solar gravity ( I varies ). 

This section is concerned with the long period (30 years) 
changes in e and I . Since these two motions are indep- 
endent (so long as e and I are relatively small) they 
will be considered separately for several types of orbits. 

All numerical results discussed in this section were ob- 
tained with the STEPR multi revolution integration program 
(References 18 and 19). Gravitational perturbation included 
geopotential terms up to order and degree of 6. The sun and 
moon were treated as point masses. The model for solar rad- 
iation pressure is described in Section 2.3. 


4.1 Eccentricity Variations 

The variations in e have been thoroughly discussed in 
Section 3. However, as mentioned there, that analysis did 
not take into account the coupling b«*'''ween solar radiation 
pressure and gravity. This effect produces a long term 
variation in e 

\ 

The accuracy oi equation (3.69) compared to a numeri- 
cally integrated solution containing gravity is shown in 
Table IV. The orbits compared are the same as in Table II. 



e and g. are the values obtained from the analyti- 
cal solution. e^ and g^ are the numerical solutions 
obtained from the STEPR program. This is a comparison 
of the force model used in Section 3, not of the analyti- 
cal solution method. 


TABLE IV. -ANALYTICAL SOLUTION VERSUS STEPR SOLUTION 



(a) 


0# 

undefined 



Years 

®A 

®S 

% 


Percent 

^A 

Error 

«A 

0.5 

.0418 

.0431 

97.3® 

99.0® 

3.0 

1.6 

1.4 

.0391 

.0392 

86.0 

9c. 9 

0.3 

5.3 

2.6 

.0405 

.0435 

129.2 

137.2 

7.9 

5.8 

4.6 

.0422 

.0460 

136.9 

151.0 

8.2 

9.3 

9.6 

.0481 

.0438 

132.^ 

165.9 

9.8 

20.0 



(b) 

®o = 

.0210 

0 
00 

1 

H 

0 

.6® 







Percent 

Error 








Years 

®A 

®S 

Sa 


®a 

«A 

0.5 

.0209 

.0222 

95.2® 

95.7® 

6.1 

0.5 

1.4 

.0193 

.0202 

71.6 

73.1 

4.6 

2.1 

2.6 

.0247 

.0250 

154.0 

154.1 

1.4 

0.1 

4.6 

.0286 

.0278 

163.3 

166.6 

2.9 

2.0 

9.6 

.0328 

.0243 

153.7 

168.5 

35.0 

8.9 
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Comments on Table IV 

1. The analytical solution is accurate to within a few per- 
cent within the first few years. Hierefore, the figures 
shown in Section 3 will be accurate to within a few per- 
cent and describe the general character of the motion. 

2. Long period effects due to gravity become loqportant 
after about eight or nine years. 

3. 'nie analytical solution could be useful for orbit pre- 
diction and control over a few years. 

Long term variations in e have been studied for the 
cases shown in Table V. Figures 10 through 13 show e versus 
time for 30 years. For each case, the inclination was approx- 
imately 7.3^. However, the motion of e was essentially the 
same as for the case when the inclination was initially zero. 


TABLE V. -INITIAL CONDITIONS FOR EXAMPLES 


Case No. 

®o 

®o 

Initial Epoch 

1 

.0 

undefined 

noon , 1 Jan . 

1980 

2 

.0210 

-80.6® 

noon, 1 Jan. 

1980 

3 

.0214 

-67.4 

noon, 1 Jan. 

1980 

4 

.0 

undefined 

noon, 3 April 

1980 


The plots shown in Figures 10 through 13 were output from 
the STEPR program, using the CALCOMP plotter hardware and 
software that is available on the Johnson Space Center Univac 
1110 con^uter system. 



Dlscuaslop of Results 


Case 1 (Figure 10) 

(1) The curve doesn't being exactly at e»0 since there 
is no output from STEPR until after 44 revolutions. 
This is because the extrapolation table needs to be 
built before the first mult irevolut ion step is taken. 

(2) The long term gravitational effects can be seen by 
conq[)aring Figure 10 with Figure 8. (Note that the 
scales are different). Both curves have the same 
general shape. The linear trend of Figure 8 has 
been moderated in Figure 10. 

(3) Large values of e can occur for an uncorrected 
SPS orbit. Such large values are probably unac- 
ceptable to the spacecraft and ground systems. 


Case 2 (Figure 11) 


This case can be compared with Figure 9. The al- 
BK>st linear incre at 'a amplitude levels off after 
about 20 years, uowever, the curves are very simi- 
lar for the first ten years. 


Case 3 (Figure 12) 

This cases corresponds to Case 4 (Figure 5(b)) in 
- - Section 3 where the value of e is limited during 

- the first few years. Also, if Figure 12 is shifted 
about an inch to the left along the time axis, it 
is almost identical to Figure 11. This can be seen 
by overlaying the two figures. Initial conditions 
of the orbit can be manipulated so that the nearly 
constant e phase occurs anywhere along the time 
axis. This may be done to achieve certain desirable 
results for station keeping purposes. It must be 
remembered, however, that since the problem depends 
on time, the epoch of initialization must also be 
taken into account. (Compare Case, 1, 5 and 6 in 
Section 3). 



Bcc«otricity 


> I 



Time (years) 

FIGURE 10: Long. pe.rlod Varlaeions In Eccentricity, Case I 
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Case 4 (Figure 13) 

The iuitial conditions are the same as for Case 1 
except that the epoch is three months later. (Com- 
pare Case 1 and Case 5 in Section 3.) The linear 
component in e is evident ov^^r the first ten years. 
However, the long period effects eventually cause a 
decrease in eccentricity so that the maximum e 
for this case is less than the maximum e for Case 
1. This example shows the importance of the long 
period gravitational terms. 


4.2 Inclination V ariations 

The motion of the orbital plane has been thoroughly dis- 
cussed in Reference 7. The discussion given here presents a 
numerically integrated sol» ion and some plots of inclination 
as a function of time. Table VI shows the cases that are dis- 
cussed here. 


TABLE VI . -INITIAL CONDITIONS FOR EXAMPLES 


Case No. 


"o 

Initial Epoch 

1 

0® 

undefined 

noon, 1 Jan. 1980 

2 

7.3® 

0.® 

noon, 1 Jan. 1980 

3 

2.0 

270.0 

noon, 1 Jan. 1980 


For each case, the eccentricity was small and had nc appre- 
ciable effect on inclination or node. 



Inclination Vusg*) 




laclinatioQ (deg.) 



Inclination (deg.) 
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Case 1 (Figure 14) 

If initially zero , the inclination increases to 
about 15* after 26 years For the first few years, 
the increase is almost linear at the rate of .859 
degrees per year. 

Case 2 (Figure 15) 

The inclination is nearly constant. The long period 
oscillation observed in this figure is caused by 
the precession of the moon's orbital plane on the 
ecliptic. This motion has a period of 18.6 years 
and depends on the epoch of initialization. 

Case 3 (Figure 16) 

When the node is initially near 270 degrees, the 
inclination will decrease to almost zero and then 
increase. With this approach, inclination can be 
kept small over a longer period of time. The ad- 
vantage is that out of plane station keeping ma- 
neuvers will be reduced or eliminated. The effec- 
tiveness of this procedure depends to some extent 
on the orientation of the moon's orbit and there- 
fore the epoch. 


5.0 Daily Effects Due to Nonzero Eccentricity a n d Inclination 

The changes in eccentricity and inclination that were 
described in the previous sections will cause a daily motion 
of the satellite, as observed from the rotating earth. It 
will wander north and south of the equator as well as east 
and west of the mean longitude. The ground track may be a 
circle, ellipse, figure eight or some other shape depending 
on the values of e and I . Usually there are restrictions 
on the allowed latitude and longitude deviations of a geo- 
synchronous satellite, because of requirements on satellite 
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and ground systems. This section will show some typical ground 
tracks for the case of small values of eccentricity and incli- 
nation. 


5.1 Development of Equations 

The frame of reference will be a rectangular coordinate 
system with x- and y-axes in the equatorial plane and 
z-axis toward the north pole. The x-axis Is directed toward 
the vernal equinox. Figures 17 and 18 define the symbols to 
be used in this development. ^ is the angle between the 

w 

Greenwich meridian and the x-axis . 

Equations will be developed here that give latitude 
and longitude X as a function of time. Begin with the ex- 
pressions for rectangular coordinates in terms of elements, 

X = r (oos n eo8 v - sin G sin v aos I) , 

y = r (sin n sin v + aos fl sin v cos I) , (5.1) 

z = r sin v sin I , 

where v = ai + f (5.2) 

z 

Since — = sin i|; 

r 

sin ip = sin v sin I (5.3) 

Note that is uniquely determined by (5.3). 

Development of equations for the longitude A are some- 
what more difficult. First, notice that expressions must be 
f nd for both sin A and ons A . From Figure 17, 

X y 

— = cos (p — = sin 4> , (5.4) 

P P 
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80 that 


sin — (ain Q oos v *** aoa Q ein v oo8 I) 

P (5.5) 

X 

ooe (|> a — (o<78 n 808 V - Bin Q sin V 008 I) 
p 


The longitude can be expressed in terms of <|> and by 


X ■ (|) - (j>^ , 


so that 


$in X “ sin 0 eo8 “ ®08 ^ sin (|> 


008 X » 008 4> 008 4>^ + sin ^ ein 


(5.6) 


Insert (5.5) into (5.6), collect texmis and make use of trigo- 
nometric identities. After introducing tbe small parameter 


a ® — (1 - cos I) , 
2 


(5.7) 


the results are 

sin X = - |(1 - a) sin (v + fl - (>^) - a sin (v - + ^^)j 

^ (5.8) 


008 


X » — " a) cos (v + ft - 4>^) + a oos (v - ft + 


where 


, , -i 

— ® (1 - sin^ I sin^ v) 


r 

P 


(5.9) 


Th^' true anomaly f is given as a function of time by 
the implicit equation 


nt - t “ E - e sin E 
o 


(5.10) 
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and 


t - 2 tan 



■ '*"> 
+ e 


1 - e 



( 5 . 11 ) 


The keplerlan elements n (mean motion), e, I, u, 9, can be 
assumed constant over a few days. 

It is desirable to measure the longitude relative to 
some "mean" value. Define the new variables 


and 


6 = X - e 


B - l + fl - <I)^^ + 0) 
o eo 


The rotation of the earth is expressed as 


A * w t + d) 

e eo 


For geosynchronous satellites, 


so that 


( 1 ) “ n 

e 


A = nt + d) 


Making use of the notation 


( 5 . 11 ) 

( 5 . 12 ) 


we have 


il = nt + 


A = n + d> — a 

®o o 


(5.13) 


Inserting (5.13) into (5.8) and using (5.12), 


in \ — Hl-a) sin (f - £. + 0) - a sin ( f + £ + 2w 

P *- 

= — 1^1-a) aos (f - £+6) + ol aos ( f + £ + 2u 

0 L 


008 \ 


e)^ , 

(5.14) 

e,] . 
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Since |6| will be less than 90®, it can be defined by 
ein 6 . Make use of 

ein 6 ■ sin X oos 0 - oos X sin 0 . (5.15) 

Substitute (5.14) into (5.15) and collect terms. The result 
is 


sin 6 


- fd - a 

n I— 


) sin (f - t) - 


a sin (f + t + 2(d) 


>] 


(5.16) 


Finally, equations (5.3), (5.9), (5.10), (5.11) and (5.16) 
define the ground track as a function of time. 


5.2 Small Eccentricity and Inclination 

It is desirable to express latitude and longitude 6 
as explicit functions of time, l.e. as functions of t 
This can be done for small e and I by making use of power 
series expansions. Assume that 

e < .064 , I< 7.3® 


Define 


1 

b = — sin I 
2 


then 


b = e 


(5.17) 


Power series expansions will be carried out in terms of e 
and b , keeping terms of order e®, b^ and eb 
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5.2.1 Longitude Equation 

From equation (5.16), Bin 5 is expressed in terms of 
f and I . Consider 


!.[x- 


I atn* (f + w) 


>] 


-i 


This expression can be rewritten in terms of b , 


:.[a- 

n t- 


2b*) + 2 b* ooe 2 (f + w) 


>] 


-i 


Expanding the above expression with the aid of the binomial 
theorem, 


r 

- * (1 - 2 b*) - b* 008 2 (f + w) + 0(b‘') 

P 

Also, 

1 

a = — (1 - ooe I) * b^ + 0(1^**) 

2 


Therefore, after truncating terms of order b** , the expres- 
sion for sin 6 is 


ijin 6 = (1 - 3b^) sin (f - t) 

1 

- “ b^ sin ( f + )?. + 2(A)) 

2 

1 

- - b^ sin (3f - i + 2(A)) 
2 
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Notice that 3b‘ can be dropped from the first term 
since Bin (f-&) ■ 0 (®)« Also, f can be replaced by I 
in the last two terms, for similar reasons. Then, 

ein 6 ■ ein (f - A) - b^ sin 2 (t + w). (5.18) 


Making use of the Equation of the Center » 

5 

f - Jl = 2e sin % + — e^ ein 21 + • • • , (5.19) 

16 


and the expansion for sine, 

1 

ein 0 = 0 03 + ... , (5.20) 

6 


the explicit expression for longitude is 

5 

6 = 2e ein £. + — e® sin 2l 
16 


- b^ ein 2 ()!. + w) + 


(5,21) 


5.2.2 Latitude Equation 

The sine of latitude is 

sin i|> = 2b sin (f + w) 

Using (5.19), 

sin (f + w) = sin (S. + w + x) 


where 

5 

X (e, H) “ 2e ein i + — e* ein 21 + •• • 

16 


.(5.22) 
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Uslng the power series expansions for sine and cosine, and 
the binomial theorem, 

008 i|) “ 1 - e* + e^ 008 2 A + * * • , 

5 

ein 1 } = 2e sin t + — e* sin 2t + • ♦ • 

16 

Then sin ip can be written as 

sin \J» = 2b sin (t + w) 

+ 2b efsin (21 + w) - sin («] + ••• 


The latitude is then 


^ - 2b ein (S. + w) 

(5.23) 

+ 2b e|sin (21 + w) - sin + Q(®**^*) 


5.3 Ground Tiack 

To an error of less than one half degrees, the second 
degree terms can be neglected in equations (5.21) and (5.23), 

6 2e sin 1 , (5.24) 

(|» = 2b sin (1 + «) (5.25) 

An approximate equation for the ground track can be obtained 
by eliminating 1 from the above equations. From (5.25) 

Ip 

— - sin 1 008 0 ) = 008 1 sin u 
2b 
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Squarlng both sides, 

— - — sin % QOS w + sin^ % jos‘ o) = aos^ £, sin^ w 
4b^ b 


But 

6 

sin t = — 
2e 


so that 

lj/^ $ t(> QOS (A) 6^ 

- + = sin^ w (5.26) 

4b^ 2b e 4e^ 


Equation (5.26) is an ellipse, so that, in general, the 
ground track will be nearly an ellipse. Consider the sp'^ial 
cases : 

( 1 ) 0 ) = 90 ** 

Then the ground track equation becomes 

4e^ 4b^ 

This is an ellipse whose axes lie- on the equator and a 
meridian . 

(2) e = b , u) = 90° 

The ground track is a circle of radius e. 

(3) 0) » 0 

The ground track is a straight line with slope 

ijj b 


6 e 
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Actual ground tracks are shown in Figures 19(a), 19(b), 
20(a), 20(b) for different values of e and I . These were 
produced on the Hewlett-Packard 9810 by using equations (5.3), 
(5.9), (5.11) and (5.16). Thus, there are no approximations 
made for small e and I . It is observed that the figures 
resemble ellipses. 


6.0 Conclusions 

The analysis developed in this report shows that the or- 
bital eccentricity of the SPS can get relatively large. How- 
ever, for certain cases, the eccentricity can be reduced when 
proper choices on iritial conditions are made. 

An analyt: solution for the motion of eccentricity 

and longitude :>j. perigee has been derived. This solution is 
valid for eight to ten years. It could be used for prediction 
and control of the SPS orbit. 

In orocr for the analytical solution to be valid for 
longer periods of cime, the gravitational effects must be 
included. It has been shown by numerical integrations that 
gravitational perturbations on the eccentricity become impor- 
tant for time intervals longer than ten years. The complete 
analytical solution is feasible through use of the methods 
developed in this report. 



Latitude (deg. 




Longitude (deg.) 
FIGURE 19; Concluded 
(b) u)*-45 , 1 35®, 18 C 
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Longitude (deg.) 

FIGURE 20: Ground Track, e*.021 



Lon gi rude ( deg . ) 
FIGURE 20: Concluded 
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APPENDIX A 

POSITION OF THE SUN IN EQUATORIAL COORDINATES 


It will be assumed that the sun moves on an elliptical 
orbit about the earth. Using the coordinate system described 
in Section 3.1, the direction cosines of the sun are: 


-2 = (f^+g@) 





C sin 


(Al) 


~ S sin > 

where C = oos e , S - sin e , 

and e is the angle between the equatorial plane and the 
ecliptic plane (€= 23*'27'). and g^ are the true anomaly 

and argument of perigee, respectively, of the sun. 

Equation (Al) needs to be expressed in terms of the sun's 
mean anomaly. This will be done by using power series expan- 
sions in the eccentricity e® of the sun’s orbit. Since e^ 
is small (eg.= .0167), only first degree terms are needed: 


f® = - + oos + e^ 2i^ 

sin fg, = sin sin 21^ 


(A2) 


Using (A2) in (Al): 
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• ooe ooe s® * «<j> ooe (2i^+g<j>) 

J*o 

y_ 

— • C ain (A^+g<j,) “ C ein * C sin (2A^+g^ 

*© 

— ^ S ain (Ae,+g©) - -5 atn g© + 5 ein (2l^^g^) 

Terms that are periodic in with coefficient e® will 

not be significant and can be neglected. The final expres- 
sions are then 

— * 008 ( ■* % ' 

» C s^n (S,^+Bc>) " c e^ ein , (A4) 

^ S ain (^<s,+g@) “ 5 e^ ain 
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APPENDIX B 

NON-SINGULAR DIFFERENTIAL EQUATIONS 


The differential equations (3.47) for p and q will 
be derived from the chain rule (3.46) and the expression 
(3.40) of Fg as a function of the DSu variables. The fol 
lowing derivatives are needed: 


3F 3 


8 _ 


dG 2 


rn- 


* — e 1 — (ooa h gob g - sin h coe I sin g) 


^ e 


- Bin h QOS I e Bin ^('?<?e v - e^ gob g^) 

r 

+ \c — iain h gob g + gob h gob I ain g) 
L e 


(Bl) 


+ S — ein I sin g 
e 

GOB^I -I 

+ C GOB h GOB I e sin g + S — ; e sin gj (gob v 

Bin I 


- e© g_) 


3F 3 ( 

— 2 - = — e j Bin h 


dH 2 


1 


e Bin g {gob V - 6gj gob 


- C 


GOB I “I 

e sin g 

Bin I J 


C coe h e sin g - S 


(B2) 

{sin \» - 


e© e®) 
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^^8 ^ r “! 

® — G £ Jle sin g ooa h + e qob g sin h cos I (cos v 

3g 2 L -I 


- cos gg,) + 

+ (e cos g sin h - e cos g cos h cos I) - 

- 5 e cos g sin ll (sin v - sin g^) 



(B3) 


Inserting these derivatives into (3.46) and collecting terms, 
the following expressions are obtained 


dT 



008 h sin h (1-coa I) 


+ (cos I-l) e^ 


sin g sin (g+h) sin h 


(cos V - 


- cos g,^) 

- [c (l+cos h (cos I-l)) + 5 


(B4) 

sin I cos h - 


- sin (g+h) (cos I-1)(C cos h sin g - 
cos I -| 

- S ) (sin V - e^ sin g^) . 

sin I -I 



-101- 


dq 

dx 


- e J (1 - (1-008 I) sin^ h) - 

2 |L 

- e^ooe (g+h)(<?o8 I-l) s^n h sin ^(oos v 

_j 

* ^ ^ (1-<20S I) 

- T)" 5 Bin I sin h 

^ cos (g+h)((?<?s I-1)(C 008 h sin g - 
cos I -] 

- S sin g; (srn v - e sin 

sin I J 


T‘ above equations can be expressed in terms of p,q 
r.^i.xng use of the following expressions; 


p = e oos (g+h) 


q = e sin (g+h) 


P = /1-oosl’ oos h 
P 


Q = - /l-<?osl’ sin h 

q 


cos ( g+h ) = 




2^ 2 


sin (g+h) = 


p ■^q 


P^+q^ 


oos h = 


/P^+(^ 

sin I cos h = P ^2 - (P^+Q^)J^ 
sin I st^ h = - Q 1^2 - (P^+Q^)J^ 
oos g = (pP-qQ) |^(P'+q^) (P^+Q^) 
sin g = (qP+pQ) |^(p^+q^) (P^+0^) 


sin h = 


-Q 


/P^+Q^’ 


-i 


e© g^) 
(B5) 


P,Q by 


(B6) 


The resulting expressions are given in (3.47). 
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The chain rule applied to the derivatives of P and Q 
result in 


4 


i 


dP eoa h 


dT 2G/1-0O8 T 


[ SF df 

- 008 I — 

3h 


3P 

sin h/i-eo8 I* — 
3H 


ir* -8in 

^ ^ m 

dT 2G/r 


in h pF 
—008 t* [_3h 


(B7) 


3F 

eo« I — 


3F 

008 h/i-eo8 l‘ — 
3H 


The following ac*ditional partial derivative is needed: 


ap 3 

— i = - e G 
3h 2 


je 008 g sin h + e sin g cos I cos ^(oos v 


6gj oos gQ> + C |e sin g oos I sin h - (BB) 

e oos g 008 b^(sin v - e^ sin g^) 


Inserting (B2), (B3) and (B8) into (B7): 


dp 

dx 


3 

- e 

4 


/1-C08 I* 1+008 I) OOS h e oos g - 

- 2 sin h e sin ^ ** (oos \> - oos 

- C oos h (sin v - e^, sin g®)] * 

+ S 008 I j^os h /l+oos r e oos g - 
- e sin g (sin v - 

r J 


2 sin h 


stn 


/l+oos 
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dQ 

dT 


3 

— e 

4 


/r 


-C 08 |( 1+008 I) sin b e cos g 


(B9) 


•¥ 2 008 h e ein 


i-n ^ 


* 008 h (ein v - e^, ein g^) - ein h (eoe v 

- ©0 fe)] - 5 e<?e I jein b /l+<?oe F e eoe g 
2 <j<3e b 1 

+ — e ein g (ein v - ^ ei« g^) 

/l+<?oe f 


The above equations can now be expressed in tereins of p.Qi 
P»Q by the use of (B6). The resulting non-singular ex- 
pressions are given in (3.48). 
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APPENDIX C - COMPUTATIONAL ALGORITHM FOR THE ANALYTICAL SOLUTION 
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APPENDIX C 

COMPUTATIONAL ALGORITHM FOR THE ANALYTICAL SOLUTION 

The approximate analytical solution for p and q was 
derived in Section 3.7 and is written as 

p « d>g (eo« 0 + T 6 e^ Bin g ) + C, 

(Cl) 

q » (sin © - T 6 e gob g ) + C„ 

® 0 4 

The computational sequence for evaluating p and q at any 
time is shown below. 

(1) Values of constants 

5 e^ sin g^ = (365.25)"^ (0.01675) sin (281®. 0) 

6 e^ 008 g^ » (365.25)“* (0.01675) OCB (281® .0) 

8 = 0.9679 , 6 = (365.25)“* 

(2) Value of small parameter 

S ^ 

e » (6.611)* - (5.06-10“’^ 

M 

S * cross-sectional surface area in square meters 

M weight in kilograms 
3 

- (365.25) e 
2 
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es) Compute Cj and C 2 

C, . e„ 00 , 2^ - $6 008 + K,) 

Cj - 0 in 2^ - $.i« (i^ + Eg) 

= initial value of eccentricity 

= initial longitude of perigee 

ft « 358® 28 ’33" + 1295 96579” T 
00 

g » 281® 13* 15" + 6189" T 
0 

T = Julian centuries of 36525 ephemeris days , 
referenced to 1900 January 0.5 , 

(4) Cong>ute p(t) and q(x) using equation (Cl), where 

T = ( 2 ir) (number of revolutions) 

e - (365.25)-* T + + gg 

It can be assumed that one revolution is equivalent to one day. 



